package fem2.element.crack_sda;

import inf.jlinalg.Array1DMatrix;
import inf.jlinalg.Array2DMatrix;
import inf.jlinalg.ArrayVector;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVector;

import java.util.ArrayList;

import math2.LinAlgUtilities;
import fem2.Debugger;
import fem2.MaterialModel;
import fem2.MeshPart;
import fem2.Node;
import fem2.material.cohesive_law.CohesiveLaw;

public class KOSConstantJumpTetrahedraElement extends SOSConstantJumpTetrahedraElement {

	public KOSConstantJumpTetrahedraElement(MeshPart mp, MaterialModel mm, CohesiveLaw cl) {
		super(mp, mm, cl);
	}

	// @Override
	// protected void computeEnhancedStrainOperator(MeshPart mp, IMatrix g,
	// double... xi) {
	// int dim = mp.getDimension();
	// int nn = mp.countModes();
	// int ntrs = dim * (dim + 1) / 2;
	// Array1DMatrix B = Array1DMatrix.createColumnMajor(ntrs, dim * nn);
	// computeB(mp, B, xi);
	// int shift = dim * posNode;
	// g.set(0, 0, -B.get(0, shift));
	// g.set(1, 1, -B.get(1, 1 + shift));
	// g.set(2, 2, -B.get(2, 2 + shift));
	// g.set(3, 0, -B.get(3, shift));
	// g.set(3, 1, -B.get(3, 1 + shift));
	// g.set(4, 1, -B.get(4, 1 + shift));
	// g.set(4, 2, -B.get(4, 2 + shift));
	// g.set(5, 0, -B.get(5, shift));
	// g.set(5, 2, -B.get(5, 2 + shift));
	// }

	@Override
	protected void computeEnhancedStrainOperator(MeshPart mp, IMatrix g, double... xi) {
		double coeff = -crackMeasure / elementMeasure;
		IVector normal = bestMode.getSeparationNormal();
		g.set(0, 0, coeff * normal.get(0));
		g.set(1, 1, coeff * normal.get(1));
		g.set(2, 2, coeff * normal.get(2));

		g.set(3, 0, coeff * normal.get(1));
		g.set(3, 1, coeff * normal.get(0));

		g.set(4, 1, coeff * normal.get(2));
		g.set(4, 2, coeff * normal.get(1));

		g.set(5, 0, coeff * normal.get(2));
		g.set(5, 2, coeff * normal.get(0));
	}

	// @Override
	// protected boolean checkOnsetCriteria() {
	// IVector sigma = getStress(0, stressType);
	// MeshPart mp = getMeshPart();
	// int n = mp.getBasis().getDimension();
	// int dim = mp.getDimension();
	// int ntrs = dim * (dim + 1) / 2;
	// /*
	// * perform spectral decomposition
	// */
	// IVector ps = new ArrayVector(dim);
	// IMatrix pd = new Array2DMatrix(dim, dim);
	// LinAlgUtilities.decomposePrincipalStress(dim, sigma, ps, pd);
	//
	// if (ps.get(2) < getCohesiveLaw().getTensileStrength()) {
	// return false;
	// } else {
	//
	// /*
	// * realize crack normal
	// */
	// double n1 = pd.get(0, 2);
	// double n2 = pd.get(1, 2);
	// double n3 = pd.get(2, 2);
	//
	// /*
	// * compute B operator
	// */
	// Array1DMatrix B = Array1DMatrix.createColumnMajor(ntrs, dim * n);
	// double[] xi = getGaussPoint(0);
	// computeB(mp, B, xi);
	//
	// /*
	// * search the gradient which make smallest angle with crack normal
	// */
	// double cos = 0.0;
	// for (int i = 0; i < 4; i++) {
	// double gradx = B.get(0, i * dim);
	// double grady = B.get(1, 1 + i * dim);
	// double gradz = B.get(2, 2 + i * dim);
	// double normgrad = Math.sqrt(gradx * gradx + grady * grady + gradz *
	// gradz);
	// gradx /= normgrad;
	// grady /= normgrad;
	// gradz /= normgrad;
	// double tmp = Math.abs(n1 * gradx + n2 * grady + n3 * gradz);
	// if (tmp > cos) {
	// cos = tmp;
	// posNode = i;
	// normal[0] = gradx;
	// normal[1] = grady;
	// normal[2] = gradz;
	// /*
	// * choose the normal making smallest angle with gradient
	// * vector
	// */
	// double tmp1 = n1 * gradx + n2 * grady + n3 * gradz;
	// normal[0] = n1 * Math.signum(tmp1);
	// normal[1] = n2 * Math.signum(tmp1);
	// normal[2] = n3 * Math.signum(tmp1);
	// /*
	// * compute crack area
	// */
	// crackArea = elementVolume * normgrad;
	// }
	// }
	//
	// /*
	// * choose tangent to align with other principal stress
	// */
	// tangent1 = new double[] { pd.get(0, 0), pd.get(1, 0), pd.get(2, 0) };
	// tangent2 = new double[] { pd.get(0, 1), pd.get(1, 1), pd.get(2, 1) };
	//
	// /******************************
	// * compute crack boundary
	// ******************************/
	// int[] negNodes;
	// if (posNode == 0) {
	// negNodes = new int[] { 1, 2, 3 };
	// } else if (posNode == 1) {
	// negNodes = new int[] { 0, 2, 3 };
	// } else if (posNode == 2) {
	// negNodes = new int[] { 0, 1, 3 };
	// } else if (posNode == 3) {
	// negNodes = new int[] { 0, 1, 2 };
	// } else {
	// throw new Error("positive node is wrong");
	// }
	//
	// /*
	// * realize three edges intersect with crack
	// */
	// int[][] edges = new int[][] { { posNode, negNodes[0] }, { posNode,
	// negNodes[1] },
	// { posNode, negNodes[2] } };
	//
	// /*
	// * set the crack center as element center and compute crack surface
	// * equation
	// */
	// double[] cen = mp.computeCenter();
	// double a = normal[0];
	// double b = normal[1];
	// double c = normal[2];
	// double d = -(a * cen[0] + b * cen[1] + c * cen[2]);
	//
	// /*
	// * compute intersection of temporary crack surface with edges
	// */
	// Node[] nodes = mp.getNodes();
	// double[][] X = new double[edges.length][];
	// for (int i = 0; i < edges.length; i++) {
	// double[] N1 = nodes[edges[i][0]].getPosition();
	// double[] N2 = nodes[edges[i][1]].getPosition();
	// double u = N2[0] - N1[0];
	// double v = N2[1] - N1[1];
	// double w = N2[2] - N1[2];
	// X[i] = LinAlgUtilities
	// .computeIntersection(a, b, c, d, N1[0], N1[1], N1[2], u, v, w);
	// }
	//
	// /*
	// * compute temporary area
	// */
	// double[] cross = LinAlgUtilities.cross(X[1][0] - X[0][0], X[1][1] -
	// X[0][1], X[1][2]
	// - X[0][2], X[2][0] - X[0][0], X[2][1] - X[0][1], X[2][2] - X[0][2]);
	// double tmpArea = 0.5 * LinAlg.norm2(dim, cross);
	//
	// /*
	// * compute scaling factor & actual crack center
	// */
	// double scale = Math.sqrt(crackArea / tmpArea);
	// double[] crackCenter = new double[dim];
	// double[] posPos = nodes[posNode].getPosition();
	// crackCenter[0] = (cen[0] - posPos[0]) * scale + posPos[0];
	// crackCenter[1] = (cen[1] - posPos[1]) * scale + posPos[1];
	// crackCenter[2] = (cen[2] - posPos[2]) * scale + posPos[2];
	//
	// /*
	// * compute crack boundary
	// */
	// d = -(a * crackCenter[0] + b * crackCenter[1] + c * crackCenter[2]);
	// for (int i = 0; i < edges.length; i++) {
	// double[] N1 = nodes[edges[i][0]].getPosition();
	// double[] N2 = nodes[edges[i][1]].getPosition();
	// double u = N2[0] - N1[0];
	// double v = N2[1] - N1[1];
	// double w = N2[2] - N1[2];
	// X[i] = LinAlgUtilities
	// .computeIntersection(a, b, c, d, N1[0], N1[1], N1[2], u, v, w);
	// crackBoundary.add(X[i]);
	// }
	//
	// if (getEchoLevel().contain(EchoLevelType.AUX4)) {
	// Debugger.watch("sigma = ", sigma);
	// Debugger.watch("crack area = ", crackArea);
	// Debugger.watch("crack normal = ", normal);
	// Debugger.watch("element volume = ", elementVolume);
	//
	// // cross = LinAlgUtilities.cross(X[1][0] - X[0][0], X[1][1] -
	// // X[0][1], X[1][2]
	// // - X[0][2], X[2][0] - X[0][0], X[2][1] - X[0][1], X[2][2] -
	// // X[0][2]);
	// // tmpArea = 0.5 * LinAlg.norm2(dim, cross);
	// // Debugger.watch("crack area 2 = ", tmpArea);
	//
	// }
	//
	// return true;
	// }
	// }

	@Override
	protected boolean checkOnsetCriteria() {
		IVector sigma = getStress(0, stressType);
		MeshPart mp = getMeshPart();
		int n = mp.getBasis().getDimension();
		int dim = mp.getDimension();
		int ntrs = dim * (dim + 1) / 2;
		/*
		 * perform spectral decomposition
		 */
		IVector ps = new ArrayVector(dim);
		IMatrix pd = new Array2DMatrix(dim, dim);
		LinAlgUtilities.decomposePrincipalStress(dim, sigma, ps, pd);

		if (ps.get(2) < getCohesiveLaw().getTensileStrength()) {
			return false;
		} else {

			alpha = new ArrayVector(dim);

			/*
			 * compute B operator
			 */
			IMatrix B = Array1DMatrix.createColumnMajor(ntrs, dim * n);
			double[] xi = getIntegrationPoint(0);
			computeB(mp, B, xi);

			/*
			 * check and filter all valid separation modes
			 */
			ArrayList<SeparationMode> validModes = new ArrayList<SeparationMode>();
			// int[][] configurations = new int[][] { { 0 }, { 1 }, { 2 }, { 3
			// }, { 0, 1 }, { 0, 2 },
			// { 0, 3 }, { 1, 2 }, { 1, 3 }, { 2, 3 } };
			int[][] configurations = new int[][] { { 0 }, { 1 }, { 2 }, { 3 }, { 0, 1 }, { 0, 2 },
					{ 0, 3 } };

			for (int i = 0; i < configurations.length; i++) {
				double gradx = 0, grady = 0, gradz = 0;
				for (int j = 0; j < configurations[i].length; j++) {
					gradx += B.get(0, configurations[i][j] * dim);
					grady += B.get(1, 1 + configurations[i][j] * dim);
					gradz += B.get(2, 2 + configurations[i][j] * dim);
				}
				double normgrad = Math.sqrt(gradx * gradx + grady * grady + gradz * gradz);
				gradx /= normgrad;
				grady /= normgrad;
				gradz /= normgrad;
				/*
				 * compute crack area
				 */
				double A = elementMeasure * normgrad;
				/*
				 * check the separation mode
				 */
				SeparationMode mode = new SeparationMode(mp, new ArrayVector(gradx, grady, gradz),
						A, configurations[i]);
				if (mode.validate()) { // by right, it must always be true
					validModes.add(mode);
				}
			}

			/*
			 * realize crack normal
			 */
			double n1 = pd.get(0, 2);
			double n2 = pd.get(1, 2);
			double n3 = pd.get(2, 2);

			/*
			 * find all in valid separation modes, the mode make the smallest
			 * angle with crack normal vector from Rankine criterion
			 */
			double cos = 0.0;
			for (SeparationMode mode : validModes) {
				IVector N = mode.getSeparationNormal();
				double tmp = N.get(0) * n1 + N.get(1) * n2 + N.get(2) * n3;

				/*
				 * this code below will create a temporary mode and check if the
				 * givem normal vector and crack area create a compatible
				 * separation mode
				 */
				// double sign = Math.signum(tmp);
				// SeparationMode tempMode = new SeparationMode(new
				// ArrayVector(n1 * sign, n2 * sign,
				// n3 * sign), mode.getCrackArea(), mode.getPositiveNodes());
				// if (Math.abs(tmp) > Math.abs(cos) && tempMode.validate()) {
				// cos = tmp;
				// bestMode = mode;
				// }

				/*
				 * this code purely accept any mode if the normal vector create
				 * a smaller angle with the mode. It may not create a compatible
				 * separation.
				 */
				if (Math.abs(tmp) > Math.abs(cos)) {
					cos = tmp;
					bestMode = mode;
				}

			}

			// assert bestMode != null : "no valid mode detected";
			if (bestMode == null) {
				Node[] nodes = mp.getNodes();
				for (int j = 0; j < nodes.length; j++) {
					Debugger.watch(nodes[j].getPosition());
				}
				// Debugger.watch("n1 = ", n1);
				// Debugger.watch("n2 = ", n2);
				// Debugger.watch("n3 = ", n3);
				// Debugger.watch("V = ", elementVolume);
				// Debugger.watch("B = ", B);
				throw new Error("no valid mode detected");
			}

			/*
			 * get the normal from the best mode to become normal of the crack,
			 * and also reverse direction if need
			 */
			IVector ns = bestMode.getSeparationNormal();
			double sign = Math.signum(ns.get(0) * n1 + ns.get(1) * n2 + ns.get(2) * n3);
			normal = new double[] { n1 * sign, n2 * sign, n3 * sign };
			crackMeasure = bestMode.getCrackArea();

			tangent1 = new double[] { pd.get(0, 0), pd.get(1, 0), pd.get(2, 0) };
			tangent2 = new double[] { pd.get(0, 1), pd.get(1, 1), pd.get(2, 1) };

			/*
			 * compute crack boundary. Firstly create a temporary separation
			 * mode
			 */
			// SeparationMode tempMode = new SeparationMode(new
			// ArrayVector(normal), crackArea,
			// bestMode.getPositiveNodes());
			// assert tempMode.validate() :
			// "the separation by actual crack must be valid";
			// // please enable assert for this line to work correctly
			// crackBoundary = tempMode.getCrackBoundary();

			SeparationMode tempMode = new SeparationMode(mp, new ArrayVector(normal), crackMeasure,
					bestMode.getPositiveNodes());
			crackBoundary = tempMode.getCrackBoundaryByScanning();

			return true;
		}
	}
}
